Five different tree species were used to select parameter sets to represent different types found around Santa Barbara.
Steps taken to select the parameters were as follows:
## summarized annually
dfyr %>%
ggplot() + geom_line(aes(x=as.factor(year), y=lai, col=code, group=code))
dfyr %>% dplyr::filter(year>2009 & year < 2020) %>%
ggplot() + geom_line(aes(x=as.factor(year), y=lai, col=code, group=code))
dfyr %>%
ggplot() + geom_line(aes(x=as.factor(year), y=plantc, col=code, group=code))
#dfday %>% ggplot() + geom_line(aes(x=yd, y=lai, col=code))
dfyr2 %>% ungroup() %>%
dplyr::filter(year ==2011 | year==2014 | year==2017) %>%
ggplot() + geom_boxplot(aes(x=code, y=lai, fill=as.factor(year))) + ggtitle("annual avg LAI")
dfyr2 %>% ungroup() %>%
dplyr::filter(year ==2011 | year==2014 | year==2017) %>%
ggplot() + geom_boxplot(aes(x=code, y=plantc, fill=as.factor(year))) + ggtitle("annual avg plantc")
# load in mikes data
rs_dat <- read.csv("../../data/alonzo_data_df.csv")
rs_dat <- dplyr::filter(rs_dat, fractional_cov > 0.7)
trees_rs <- rs_dat %>% dplyr::filter(
class==81 | #QUAG
class==74 | #PLRA
class==73 | #PIUN
class==66 | #PICA
class==34 ) #EUGL
## remove outliers
trees_rs_rm <- trees_rs[!trees_rs %in% boxplot.stats(trees_rs$carbon)$out,]
ggplot(trees_rs_rm) + geom_boxplot(aes(x=as.factor(class), y=lai, fill=as.factor(class))) + ggtitle("LAI estimates from alonzo 2016")
## Warning: Removed 2143 rows containing non-finite values (stat_boxplot).
trees_rs_rm$code <- NA
trees_rs_rm$code[trees_rs$class==81] <- 'quag'
trees_rs_rm$code[trees_rs$class==74] <- 'plra'
trees_rs_rm$code[trees_rs$class==73] <- 'piun'
trees_rs_rm$code[trees_rs$class==66] <- 'pica'
trees_rs_rm$code[trees_rs$class==34] <- 'eugl'
ggplot() +
geom_boxplot(data=dfyr2, aes(x="model", y=lai, fill='rhessys'),position=position_dodge2(padding=0.5)) + facet_wrap('code') +
geom_boxplot(data=trees_rs_rm, aes(x="data", y=lai, fill='alonzo'), position=position_dodge2(padding = 0.5)) + facet_wrap('code') +
ggtitle("compare LAI")
## Warning: Removed 2143 rows containing non-finite values (stat_boxplot).
ggplot() +
geom_boxplot(data=dfyr2, aes(x="model", y=plantc, fill='rhessys')) + facet_wrap('code') +
geom_boxplot(data=trees_rs_rm, aes(x="data", y=carbon, fill='alonzo')) + facet_wrap('code', scales='free') +
ggtitle("compare plant carbon (kgC/m^2)")
## Warning: Removed 2143 rows containing non-finite values (stat_boxplot).
# ggplot() +
# geom_boxplot(data=dfyr2, aes(x="model", y=height, fill='rhessys')) + facet_wrap('code') +
# geom_boxplot(data=trees_rs, aes(x="data", y=height, fill='alonzo')) + facet_wrap('code') +
# ggtitle("compare height (m)")
# ggplot() +
# geom_boxplot(data=dfyr2, aes(x="model", y=rootdepth, fill='rhessys')) + facet_wrap('code')
# load in David's data
#vis <- read.csv("/Volumes/GoogleDrive/My Drive/plant_params/DM_data/tree_and_turfgrass_monroe2_selected_weightedmeanVIs.csv")
vis <- read.csv("../../data/tree_and_turfgrass_monroe2_selected_weightedmeanVIs.csv")
trees_vi <- vis %>% dplyr::filter(sp_class==81|
sp_class==74|
sp_class==73|
sp_class==66|
sp_class==34)
vi_dens <- ggplot(trees_vi) + geom_density(aes(x=NDVI, col=as.factor(year))) + facet_wrap("sp_code")
rh_dens <- dfyr2 %>% ungroup() %>%
dplyr::filter(year ==2011 | year==2014 | year==2017) %>%
ggplot() + geom_density(aes(x=lai, col=as.factor(year))) + facet_wrap("code")
grid.arrange(vi_dens, rh_dens)
model_pc <- dfyr2 %>% ungroup() %>%
dplyr::filter(year ==2011 | year==2014 | year==2017) %>%
pivot_wider(values_from=lai, names_from=year, id_cols=c(run,code)) %>%
mutate(decline = (`2014`-`2011`)/`2011`) %>%
mutate(recover = (`2017`-`2011`)/`2011`) %>%
mutate(data="rhessys LAI") %>%
dplyr::select(ID = run, sp_code=code, decline, recover, data)
model_pc$sp_code = toupper(model_pc$sp_code)
index_pc <- trees_vi %>%
dplyr::filter(year ==2011 | year==2014 | year==2017) %>%
pivot_wider(values_from=NDVI, names_from=year, id_cols=c(ID,sp_code)) %>%
mutate(decline = (`2014`-`2011`)/`2011`) %>%
mutate(recover = (`2017`-`2011`)/`2011`) %>%
mutate(data="NDVI") %>%
dplyr::select(-`2011`, -`2014`, -`2017`)
perc_changes <- rbind(model_pc, index_pc)
dec <- ggplot(perc_changes) +
geom_boxplot(aes(x=sp_code, y=decline, fill=data))
rec <- ggplot(perc_changes) +
geom_boxplot(aes(x=sp_code, y=recover, fill=data))
perc_changes %>%
pivot_longer(cols=c(decline, recover), names_to="change", values_to="perc_diff") %>%
ggplot() +
geom_hline(aes(yintercept=0), linetype="dashed") +
geom_boxplot(aes(x=change, y=perc_diff, col=data), fill='grey') +
facet_wrap('sp_code') +
theme_bw() + theme(legend.position="top") +
labs(y="percent difference to 2011 values")
Llorett et al. 2011 defines resilience as the capacity to reach pre-episode growth levels and relative resilience as the resilience weighted by the damage incurred during the episode.
By looking at time series as each year being relative the the pre-drought conditions, we can see how many years it takes to recover (or not at all). can also figure out which year to consider the drought occurring - where do we assign drought vs. post drought?
\(resistance = Drought / PreDrought\)
\(recovery = PostDrought / Drought\)
\(resilience = PostDrought / PreDrought\)
\(RelativeResilience = (PostDrought-Drought)/PreDrought\)
In the model, tree LAI responds the worst to the drought in 2015, after the months of experiencing the extreme heat in 2014 and 2015.
p_dec <- read.csv("decid_params.csv")
p_ev <- readRDS("evg_broadleaf_params.rds")
p_con <- readRDS("pica_conifer_params.rds")
p_dec2 <- cbind(p_dec[,15:54], p_dec$code)
p_ev2 <- cbind(p_ev[,14:53], p_ev$code)
p_con2 <- cbind(p_con[,14:53], p_con$code)
p_names <- c(
"pore_size_index",
"psi_air",
"run",
"epc.leaf_cn",
"epc.branch_turnover",
"epc.gl_smax",
"epc.flnr_age_mult",
"epc.litter_moist_coef",
"epc.psi_close",
"mrc.q10",
"epc.vpd_close",
"epc.leaf_turnover",
"epc.froot_turnover",
"epc.storage_transfer_prop",
"epc.height_to_stem_coef",
"epc.waring_pa",
"epc.root_distrib_parm",
"epc.proj_sla",
"epc.alloc_stemc_leafc",
"epc.ext_coef",
"specific_rain_capacity",
"epc.flnr_sunlit",
"epc.flnr_shade",
"epc.netpabs_sunlit",
"epc.netpabs_shade",
"epc.netpabs_age_mult",
"epc.cpool_mort_fract",
"epc.min_percent_leafg",
"epc.livewood_turnover",
"epc.waring_pb",
"epc.max_storage_percent",
"epc.min_leaf_carbon",
"epc.resprout_leaf_carbon",
"epc.alloc_frootc_leafc",
"epc.frootc_crootc",
"epc.alloc_prop_day_growth",
"epc.day_leafon",
"epc.day_leafoff",
"epc.ndays_expand",
"epc.ndays_litfall", "code")
colnames(p_dec2) <- p_names
colnames(p_ev2) <- p_names
colnames(p_con2) <- p_names
all_p <- rbind(p_dec2, p_ev2, p_con2)
###### now can use to plot
ggplot(all_p) + geom_boxplot(aes(x=as.factor(code), y=epc.proj_sla, fill=as.factor(code))) + ggtitle("proj_sla")
ggplot(all_p) + geom_boxplot(aes(x=as.factor(code), y=epc.leaf_turnover, fill=as.factor(code))) + ggtitle("leaf_turnover")
ggplot(all_p) + geom_boxplot(aes(x=as.factor(code), y=epc.waring_pa, fill=as.factor(code))) + ggtitle("waring_pa")
ggplot(all_p) + geom_boxplot(aes(x=as.factor(code), y=epc.waring_pb, fill=as.factor(code))) + ggtitle("waring_pb")
ggplot(all_p) + geom_boxplot(aes(x=as.factor(code), y=epc.branch_turnover, fill=as.factor(code))) + ggtitle("branch_turnover")
ggplot(all_p) + geom_boxplot(aes(x=as.factor(code), y=epc.froot_turnover, fill=as.factor(code))) + ggtitle("froot_turnover")
ggplot(all_p) + geom_boxplot(aes(x=as.factor(code), y=epc.root_distrib_parm, fill=as.factor(code))) + ggtitle("root distrib")
comb <- left_join(dfyr2, all_p, by=c('code','run')) %>%
dplyr::filter(year==2011 | year==2014 | year==2017) %>%
group_by(year, run, code) %>% summarize_all(funs(mean), na.rm=T)
## Warning: Column `code` joining character vector and factor, coercing into
## character vector
key parameters for lai and plantc:
ggplot(comb) + geom_jitter(aes(x=epc.branch_turnover, y=plantc, col=as.factor(year))) + facet_wrap('code') +
ggtitle("branch_turnover")
## Warning: Removed 288 rows containing missing values (geom_point).
ggplot(comb) + geom_jitter(aes(x=epc.leaf_turnover, y=lai, col=as.factor(year))) + facet_wrap('code') +
ggtitle("leaf_turnover") + geom_smooth(aes(x=epc.leaf_turnover, y=lai, col=as.factor(year)))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 288 rows containing non-finite values (stat_smooth).
## Warning: Removed 288 rows containing missing values (geom_point).
ggplot(comb) + geom_jitter(aes(x=epc.waring_pa, y=lai, col=as.factor(year))) + facet_wrap('code') + ggtitle("waring pa")
## Warning: Removed 288 rows containing missing values (geom_point).
ggplot(comb) + geom_jitter(aes(x=epc.proj_sla, y=lai, col=as.factor(year))) + facet_wrap('code') + ggtitle("proj_sla")
## Warning: Removed 288 rows containing missing values (geom_point).
ggplot(comb) + geom_jitter(aes(x=epc.psi_close, y=lai, col=as.factor(year))) + facet_wrap('code') +
ggtitle("psi_close")
## Warning: Removed 288 rows containing missing values (geom_point).